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ABSTRACT 

The helioseismic observations of the internal rotation profile of the Sun raise ques- 
tions about the two-dimensional (2D) nature of the transport of angular momentum 
in stars. Here we derive a convective prescription for axisymmetric (2D) stellar evo- 
lution models. We describe the small scale motions by a spectrum of unstable linear 
modes in a Boussinesq fluid. Our saturation prescription makes use of the angular 
dependence of the linear dispersion relation to estimate the anisotropy of convective 
velocities. We are then able to provide closed form expressions for the thermal and 
angular momentum fluxes with only one free parameter, the mixing length. 

We illustrate our prescription for slow rotation, to first order in the rotation rate. 
In this limit, the thermodynamical variables are spherically symetric, while the an- 
gular momentum depends both on radius and latitude. We obtain a closed set of 
equations for stellar evolution, with a self-consistent description for the transport of 
angular momentum in convective regions. We derive the linear coefficients which link 
the angular momentum flux to the rotation rate (A- effect) and its gradient (a-effect). 
We compare our results to former relevant numerical work. 
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1 INTRODUCTION 

Computations in stellar evolution have generally been car- 
ried in one-di mensional (ID ) fram eworks. Since the sem- 
inal paper of iBohm-Vitensel (|l958T) . mixing length theory 
(MLT) has proved to be a very powerful tool to compute 
the transport of heat in stars, even though its underlying 
assumptions are very often regarded as crude in comparison 
to the complexity of the usually highly turbulent convective 
flows. Recently, it has also become clear from helioseismic 
observations that the rotation al profile of the Sun is intr insi- 
cally two-dimensional (2D, see lSchou fc co authordll998l . for 
example). Moreover, the inclusion of rotation in stars has 
been s hown to be essential in many phases of stellar evo- 
lution (|Mevnet fc Maederl Il997l ; lYoon fc Langerl 120041 ') but 
these simulations usually assume solid-body rotation within 
convective regions and a self-consistent treatment of rota- 
tion and convection is still lacking. We shall emphasise in 
this work that the characteristics of the angular momentum 
fluxes depend on the latitude even for spherically symmetric 
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rotating stars. Indeed, the properties of the turbulent mo- 
tions should naturally depend on the angle between the grav- 
itational field and the angular velocity vector. It is thus clear 
that a proper treatment of the evolution of the rotation pro- 
file of stars requires a two-dimensional description. In this 
work we lay the basis of a self-consistent MLT formulation 
for 2D-axisymmetric rotating stars which could be adapted 
for future 2D stellar evolution computations. Wc illustrate 
our precriptions in the case of slow rotation, to first order 
in the rotation rate. In this limit we recover the classical 
ID stellar evolution set of equations for spherically symetric 
stars with an additional 2D equation for angular momentum 
transport. We also provide a spherically averaged version for 
practical uses in current ID stellar evolution codes includ- 
ing slow rotation, which allows for a self-consistent treat- 
ment of angular momentum transport in convective regions. 
We discuss our results in comparison to appropriate existing 
numerical simulations. 

We start from the stellar fluid dynamical equations 
which we average on a smoothing length scale. Next, we 
consider the equations for perturbed quantities, for which 
we use a quasi-linear approximation: namely, we assume the 
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perturbed fields are composed of a spectrum of unstable 
linear modes. Their amplitude is then determined by a sat- 
uration prescription which takes into accounts some of the 
non-linearities. We finally proceed to compute the heat and 
momentum fluxes which enter our averaged equations, thus 
closing our system. The only parameter in our model is the 
smoothing length scale, which we identify with the mixing 
length. This enables us to construct the stellar evolution- 
ary equations without invoking additional parameters with 
accompany ing assumptio ns. 

Earlier iGougbl \l97$F \ and lDurnev fe SpruitJ <|l979l) de- 
veloped very similar ideas but they adopted slightly differ- 
ent saturation prescriptions. They modelled the perturbed 
quantities with a single representative unstable mode 
with unspecified parameters to characterize its orientation 
whereas here we use the linear dispersion relation to infer 
the full spectrum of the perturbations. We therefore predict 
the anisotropy without extra parameters. An advantage of 
our approach is that it spells out the underlying assumptions 
which can the n form the basis to improve the prescription. 
lOrilviel l|2003h followed bv lGaraud etaH (|2010t ) derived dy- 
namical equations for the second order correlations supple- 
mented by closure relations which introduc e a set o f addi - 
tional non-dimensional parameters. Earlier, ICanutol (|l997l ) 
went even further in the hierarchy of correlations and pro- 
vided a set of dynamical equations for quantities up to third 
order correlation terms with closure relations param e trized 
by even more free parameters. iKichatinov fc Rudigerl (|l993l ) 
approximated the effects of turbulence by a viscous stress 
tensor and computed the effects of rot ation on the mo- 
mentum fluxes. Except for I Canute! l|l997l ). all these authors 
discussed only solid body rotation. We consider here the 
dynamics of convective motions in the presence of large- 
scale fields, such as non-uniform rotation and incorporate 
the presence of a shear in the calculation of turbulent fluxes. 

In Section 2 we describe our framework in a Cartesian 
geometry and we present and discuss our MLT pre scrip- 
tion, comparing i t with the works of iGoughl (|f 9781 ) and 
iDurnev fc SpruitJ (|l979h . We apply it to slowly rotating 
spherical axisymmetric stars in Section 3. We discuss our 
results and conclusions in Section 4. In the Appendix, we 
give the full closed set of stellar evolution equations in the 
limit of slow rotation and we provide their spherically aver- 
aged equivalent for ID stellar models. 



dt(pvi) + dj(pViVj +p5i 



P9i 



(2) 



where p is the pressure and we have neglected viscosity. 

The heat transp ort equation (see for example 
lLandau fc Lifshit3l 19871 . equation 49.4) combined with con- 
tinuity (equation [TJ becomes 



Tdt(pS) + Tdj{v jP S) + OjF, = q, 



(3) 



where S is the specific entropy, q the net heat generation 
rate, T the temperature and it is assumed there is no depen- 
dence of entropy on the chemical composition. The radiative 
flux is given by 



Fj = 



-pepxdjT 



(4) 



where c p is the heat capacity at constant pressure, \ i s 
the thermal diffusivity (unit length x velocity). The spatial 
derivatives of both c v and \ are assumed to be negligible. 



2.1 Average equations 

We now define the sliding average of a quantity y at position 
x by 



yd x 



(5) 



where V(x) is a small cube of volume V centred on position 
x. We want to find a new set of equations for the averaged 
quantities. These constitute our stellar model equations. We 
use as new variables the volume and mass weighted averages, 



and 



■{p), 

(pvj) 
(P) 



= (Pe)/(P)- 



(6) 



(7) 



(8) 



We then define the residuals with respect to these averages 
by 



y = y-y 



(9) 



2 GENERAL FRAMEWORK 

We start with the equations of fluid dynamics subject to a lo- 
cal gravitational acceleration g with Cartesian components 
Qi. The mass conservation is given by 



d t p + di(pVi) = 0, 



(1) 



where p is the mass density, Vi the components of the veloc- 
ity and dt = d/dt and di = djdxt are the partial derivatives 
with respect to time and each of the three spatial coordi- 
nates. We have used Einstein's summation convention. 
The momentum conservation leads to 



IGougbl i2012h extended his earlier work to higher order in O. 



for any quantity y. 

We now make use of the Boussinesq approximations 
that velocities are small compared to the sound speed, 
wavelengths are small compared to the local scale height 
and p' ~ except wh en coupled with the gra vity gi (cf. 
ISpieeel fc Veronijl96Cl ). Although lGougbl l| 19691 ) has shown 
these approximations to be not suitable for stellar convec- 
tive regions, and instead the anelastic approximation should 
be used, we feel they capture the essential physics while sim- 
plifying the derivation. In particular, these approximations 
allow us to work with nearly incompressible equations and to 
have {y') = for most quantities of interest. With these ap- 
proximations, the volume averaged continuity is unchanged, 



d t p + di(pvi) = 0. 



(10) 
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The momentum equation becomes 

dt(pvi) + dj(pViVj + pSij + U i:j ) = pgt 



(11) 



where we discard {p'g'}, because gravity is slowly varying 
and we use the convective momentum flux (more commonly 
referred to as the Reynolds stress tensor) 

Hij=p{vWj). (12) 

Finally the average entropy evolution equation is 

Td t (pS) + IVApr.S + T,) + d 3 F 3 = q (13) 

where, noting that (f + T'y 1 w T _1 (l - T'/T) we neglect 
the non-linear terms (q'T 1 ) and (T'djFj), as is necessary 
to recover the classical formulation of MLT, and define the 
convective heat flux 



Fj = piS'v'j). 



(14) 



We drop the over bars in what follows to ease writing and 
reading but they should be assumed in the remainder of this 
section. In order to get expressions for the convective fluxes, 
we now turn to the estimation of the small scales (or y ) 
quantities. 



2.2 The sub-grid model 

In order to make progress, we make assumptions about the 
linearity and scale of the perturbations. These are those usu- 
ally made for a local mixing length theory. In the linear the- 
ory, we will then apply some kind of closure condition to 
determine the amplitudes of turbulent modes. 

We appeal to the Boussinesq approximation and com- 
pute the difference between the general equations and their 
averages in order to obtain governing equations for the y' 
quantities. The continuity equation for the perturbed veloc- 
ity becomes 



divl = 0. 
The momentum equations yield 



(15) 



dtipvD+djlpviSvj+pvjv'i+pv'iv'j+SpSij-Hij] = p'gt, (16) 

where we have discarded g' because the gravitational field 
is produced by mass deep inside the star. 

We now resort to a length scale separation hypothesis: 
averaged quantities are assumed to be nearly uniform in the 
local volume V and perturbed quantities vary over scales 
that are much smaller than the diameter of V. Hence the 
gradient oilZij, 



djlZij < dj[pViv'j + pVjv'i + pv'^ +p'6i : 



(17) 



the gradients of the perturbed quantities. Mixing-length the- 
ories which make use of this approximation are generally 
called local mixing length theories. Most MLT used for prac- 
tical purposes in stellar evolution are of this type. 

We further neglect the non-linear terms and so set 



pv'iVj 







(18) 



and retain only the linear approximation. In accord with 
the Boussinesq approximation we also neglect the temporal 
and spatial variation of the average mass density. Hence, the 



perturbed quantities are determined by solving the linear 
problem for the scales that fit well inside the local volume 
V so that 



dtv'i + v'jdjVi + Vjdjv'i + -dip 



-9i 



(19) 



where we have retained the shear term from the background 
velocity. This is necessary for the redistribution of angu- 
lar momentum. With the Boussinesq approximation and ne- 
glecting the pressure perturbations with respect to thermal 
effects, we write the first law of thermodynamics as 



<?' ( T ' y? p \ T ' 

S = Cp{— - Va — ) ~ C p — 
1 P 1 



(20) 



where V a = (d logT/d log p)s is the usual adiabatic gradi- 
ent. The entropy conservation equation is then linearised as 



T' fT'\ 1 f T' 



(21) 



(22) 



where we have neglected q' , the time-dependence of c v and 
the stratification in the thermal diffusion term. As is custom- 
ary in mixing length theories, we have neglected the term 

pT\q p T 

We simply note here that this term could become important 
when strong nuclear burning takes place within convective 
regions. 

2.3 Saturation and amplitude of the linear modes 

We restore some of the non-linear effects by adopting a 
strong assumption for the saturation of each mode. We de- 
note by A m the smoothing length-scale, a typical scale of the 
smoothing volume V. At the largest scales within this vol- 
ume, i.e. scales on the order of or just below the smoothing 
length scale A m , we assume that the saturation of a given 
mode is due solely to its own shear. Parasitic instabilities, 
such as Kelvin-Helmholtz rolls, feed on the shear motions 
generated by the parent mode. We assume that eventually 
they are responsible for its saturation. We designate v\ to 
be the complex amplitude of the Fourier mode of the veloc- 
ity perturbation associated with a wave vector k and define 
the amplitude of the velocity perturbation 



u k = 



\/|w'lfc| 2 + |«'2fe| 2 +|w'3fe| 2 . 



(23) 



We write schematically the time evolution of this amplitude 
by 



0~kU k — CTp(fe, U k )u h 



(24) 



where a k is the real part of the linear growth rate of the 
mode and a p (k,Uh) is the real part of the growth rate of the 
parasitic mode responsible for its saturation and we assume 
that cr p depends only on k and on the amplitude u k of the 
parent mode. For example, in the case of Kelvin-Helmholtz 
rolls, tip depends on the component of A; orthogonal to v'k 
but, because we have assumed that the motions are almost 
incompressible, <r p is equal to kuk where k = \k\ is the wave 
number. At saturation, equation (|24p allows to write 



Ofc = <Tp(fe, Uk) 



(25) 
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which determines the mode's amplitude. For example, if 
Kelvin-Helmholtz is the dominant parasitic instability, the 
saturation is reached when the velocity is of the order of 



Uk = crk/k. 



(26) 



Note that, provided <7k depends on the direction of the wave 
vector fc, this prescription leads to an anisotropic amplitude 
of the velocity. 

Similar ideas have been used to assess the satura- 
tion of other instabilities. The parasitic instabilities of 
the magnetorota t ional insta bility (MRI) were d escribed by 
iGoodman fc Xul (|l994h and lLatter et al.l ((2009). Their role 
for the saturation of the MRI was examined independentl y 
by iLesaffre et all |2009l ) and iPessah fc Goodman! (|2009h . 
The prescription we use he re is very similar to that used by 
IPessah fc Goodman! l|2009l ) and lPessahl (|20ld ).l Guilet et al.l 
( 20101 ) have used a slight refinement of these prescriptions 
to predict the saturation amplitude of the standing accre- 
tion shock instability (SASI) in core-collapse supernovae. 
However, these ideas have so far been concerned with in- 
dividual unstable modes. We propose here to extend this 
prescription to a whole spectrum of modes. At the largest 
scale we use the prescription (|26|) . But the smallest scales 
are likely to feel the non-linear interactions of the scales im- 
mediately above an d below, as in the Kolmogorov cascade 
l|Kolmogorovl ll94lT ) . We therefore use a power-law scaling 
for each direction individually as a first approximation. We 
write fc m = 7r/A m the minimum wave number in our smooth- 
ing volume. The resulting closure expression is then 



(27) 



where fc' = with n = 11/6 for Kolmog orov scaling or 
n = 21/1 for Bolgiano-Obukhov scaling (|Bolgiand fl959l; 
IObukhovlll959l ). These will l ikely bracket the actual spec- 
trum index (see lRinconl l2006) . The last factor accounts for 
the anisotropy of the driving instability and the first for the 
energy cascade. Our closure equation completely determines 
the amplitude of all Fourier coefficients of the velociy per- 
turbations. Once the velocity amplitude is known, the linear 
system of equations for the perturbations is used to esti- 
mate the amplitude of all other perturb ed variab l es rel ative 
to the velocity. The numerical study of lRinco"nl (|2006l ) has 
carefully examined the anisotropy and scaling of turbulent 
convection and we plan to validate our approach with such 
numerical studies. Note that the fluxes involve integrals of 
u k d 3 k oc k 2 ~ 2n dk for k ^ fc m , so our results are only weakly 
sensitive to the exponent n as long as it is strictly greater 
than 3/2, otherwise these integrals diverge. This means that 
the fluxes are dominated by the largest scales just below the 
mixing length. In fact this conflicts with our assumption of 
scale separation but this is a common inconsistency of mix- 
ing length theories. 

Others have avoided such divergent behaviour by con- 
sidering only a limited number of modes. For example, in 
iGouehl 's { 1978) statistical picture eddies of a given shape 
are randomly formed, grow and get disrupted whereas we 
envisage the sub-grid scale motion s to be a collection of 
saturated unstable modes. iGoughl has an elaborate time- 
dependent model for the evolution of an eddy, whereas our 
work assumes a steady-state which saves us from specifying 



the initial conditions. For example, he has to assume seeds 
for the eddies to be isotropically distributed. Further, he 
has to prescribe a space filling factor for the shape of his 

representative eddy. 

From equation (4.6) of IGoughl (|l978l ). his definition of 
£ = n/k v and if we identify his w 2 with the square of the 
magnitude of the radial component of our velocity |i>',-fc| 2 
and his a to our ak, we arrive at 



|w'rfc| 2 = An 2 ai(kl/k 4 ), 



(28) 



where k^ is the magnitude of the horizontal part of the repre- 
sentative wave vector and A is a calibrateable dimensionless 
constant which incorporates the anisotropy parameter as 

well as the filling facto r. 

iDurnev fc Spruitl l| 19791) have also proposed a similar 
saturation prescription to ours. They set 



(«'r) 



(29) 



for a linear combination of a few modes with similar wave 
vectors. Both these prescriptions slightly differ mathemati- 
cally from ours. However, the fundamental difference lies in 
the existence of a parameter which prescr ibes the aniso tropy 
of the velocity fie l d in b oth the works of IGoughl (|l978l ) and 
IDurnev fc Spruitl (Il979l ) , whereas our prescription links this 
anisotropy to physics of the underlying instability which gen- 
erates the perturbations. 



2.4 Computation of the fluxes 

We have hitherto discussed how to determine the modulus 
of all the Fourier coefficients of the perturbations. We now 
summarize how to compute the convective fluxes, which de- 
pend on the volume average of a product of two perturbed 
quantities: 



(y 



z) = vJ v yzd 



(30) 



We assume the volume V is a cube of side A m , hence 
V = X m . Any field on this cube can be represented by its 
Fourier modes with wave number coordinates as multiples 
of fc m = 7r/A m , for instance: 



tf \ \ / ik.; 

fc/fc m ez 3 



(31) 



where y' denotes the Fourier transform of y' . We use Parse- 
val's theorem to write 



/ y'z" d 3 x = X m V V'k z '*k 
Jv , ,, „„, 



fc/fe m ez 3 



and 



/*2:'d 3 3; = A^ ^2 y'lz'k- 



(32) 



(33) 



fe/fc m ez 



We take the average of the two previous equations and use 
the fact that y' and z' are real fields to get 



Vd 3 x = A^ Y, |j/ fe ||?fe|cos(V>(fc)) (34) 



fe/fc m ez 3 



where ^(fc) is the phase difference between y' k and z' k . We 
estimate this phase difference from the linear analysis of the 
corresponding mode. 
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Finally, we approximate the sum on all wave vectors by 
a continuous integral over the non-dimensional wave- vector 



{y' z ') 



|y'fclk'fc|cos(^(fc))d 3 fc. 



(35) 



3 SLOWLY ROTATING AXISYMMETRIC 
STARS 

The local rotation rate Q of a star can be compared to 
two rates of interest to construct dimensionless numbers. 
On the one hand, the inverse of the free-fall time scale 
yields eo = Q-^/r/g, while the buoyancy frequency provides 
another number e = Q/JV where N is the magnitude of 
the Brunt- Vaisala, frequency. The latter can also be writ- 
ten e — yj Hs I g where Hs is on the order of the entropy 
scale height. Although most stars have eo <S 1 5 the entropy 
mixing in convective regions can make Hs very large and 
e is not necessarily close to zero. For example, our Sun has 
eo = 7.4 x 10 -4 at the surface, but e can be of order unity 
at the bottom of the convective region. In the following, we 
derive the stellar evolution equations to first order in the 
parameter e and assume eo < e. We will further assume that 
the tides are weak and that an axisymmetric model about 
the rotation axis can suffice. 



3.1 Background state 

In such an axisymmetric star, it is natural to use a spher- 
ical coordinate system, with r, 9 and <j> as the radius, co- 
latitude and azimuth respectively, for the definition of the 
background. The hydrostatic pressure balance with cen- 
trifugal acceleration necessarily implies that the deviation 
from spherical symetry in the thermodynamic quantities p 
and p is on the order of e^. In the first order we can there- 
fore safely assume that the thermal background depends on 
the radius r only and that the gravitational acceleration is 
radial. Similarly, since meridional circulation is the result of 
second order terms (fi 2 ), we neglect it. Then the averaged 
velocity profile consists only of cylindrical rotation so that 



and 



U0 = rQ(r, 9) sin 9 



v r = vg — 0. 



According to our assumptions the background must be 
smooth over the scale of the volume V. This requires that 
the first derivative of Q, with respect to the cylindrical radius 
is zero on the axis of symmetry. 



3.2 Linear System of equations for the modes 

We develop the perturbation at a position xq in terms of lo- 
cal Fourier modes in a local Cartesian frame rotating about 
the axis of symmetry with angular velocity Qq = Q(xo)- 
The three axes of the frame x, y and z are made to coincide 
with the local spherical coordinate unit vectors, r, and ifi 



at xo- We consider only one single mode in this subsection. 
Thus 



y h = $t[y' k exp(st + ik x x + ik y y + ik z z) 



(36) 



where y' k is the complex amplitude of the mode under con- 
sideration, i = v— 1 , x, y and z are the coordinates in this 
local Cartesian frame and k x , k y and k z are the three compo- 
nents of the wave vector in this frame. The three coordinates 
have the physical dimension of the radius r. In particular, 
dx — dr, dy — rod9 and dz — ro sin 9od<j>. In principle, shear 
deforms non-axisymmetric perturbations on a time scale of 
the order of the local shear time, which we therefore assume 
to be long compared with the growth time 5R(s) _1 in order 
to apply our linear analysis. Although this is true for the 
Sun now, it may not necessarily hold for all stars. It cer- 
tainly is not true for Keplerian discs, where the shear rate 
is comparable to the rotation rate and to the inverse of the 
vertical convective turnover time-scale. We shall drop the fc 
subscripts from the complex amplitudes in this section and 
the next in order to ease the readability. 

The linearised continuity equation (|15|l leads to the in- 
compressibility condition 



k x v' r + kyv'e + k z v'<f, = 0. (37) 

The momentum equation (|19[1 now includes an additional 
Coriolis term 2f2o X v' because the background velocity in 
the rotating frame is = r sin9[Q(r, 9) — Qo] ■ We assume 
that the apparent gravitational field (including centrifugal 
acceleration) is vertical and write g — g r , in accordance with 
our first order expansion in the rotation rate. The linearised 
Euler equations become 



2fio sin #o v' , 



%k x i 

P 

P 



? 
! — ■ 
P 



sv's — 2fio cos#o v'd> ~ -p 

P 



(38) 



(39) 



and 



where 



and 



sv',,, + 2Q. M r v' r + 2Q, M e v' e = -— p, (40) 

P 



M r 



Mo 



1 



2Qr sin 9o 



d r j 



-9e j 



2Q.r 2 sin t _ 
and the specific angular momentum j is 

j — r 2 Q. sin 2 9q. 



(41) 



(42) 



(43) 



Here, M r and Me are two dimensionless quantities propor- 
tional to the spherical coordinates of the gradient of specific 
angular momentum at our reference point so that 

1 



M = 



2Q,r sin 9, 



■Vj. 







(44) 



For uniform rotation, this vector is simply M = R where 
R is the cylindrical radius unit vector. 
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The Boussinesq approximation means that pressure 
perturbations are negligible when compared to density and 
thermal perturbations. So the equation of state becomes 



preferred. Our saturation prescription based on the direc- 
tional dependence of the growth rate will be sensitive to 
this. 

The first order of the largest real root is 



P A i- 

— + A — 

P T 



0. 



where 



A 



<91np 
dlnT 



(45) 



(46) 



si = (k.cf>)(R- M).[k + 



k.O 



1 - (k.r) 



;<>]■ 



(52) 



is the compressibility at constant pressure. 

Finally, the entropy equation (|2f I) becomes 



N 2 V' r , 



where 



(47) 



(48) 



is the square of the modulus of the wave vector and the 
thermal Brunt- Vaisala frequency TV 2 is 



N 2 = 

Cp 



(49) 



where we neglect the latitudinal thermal gradients, in ac- 
cordance with our first order expansion in the rotation fre- 
quency. In the following, we drop the subscripts of 8 and 
M for the sake of tidiness. 



3.3 Dispersion relation and growth rate 



The set of linear equations (|37[) to (|47|l forms an eigenvalue 
problem for s. Its dispersion relation is a cubic in s which 
we express as: 

s 3 + [l - (k.r) 2 ]N 2 s 

+2n(k.cf>)(R - M).[s 2 fc + (k.0)N 2 0] + 4fi 2 (fc.fi)fc.(M x 4>) s 

+Xk 2 [s 2 + 2Qk.(R - M)s + 2Q. 2 (k.£l)k.(M x </>)] = 

(50) 

where k and ft are the unit vectors along fe and tt. Without 
thermal diffusion, this dispersion relation depends only on 
the direction of the wave vector k and not on its magnitude. 
For uniform rotatio n (JVf = R) a nd y = 0, we recover 
the re sults from both lCowlind (|l95ll ) and lDurnev fc Spruitl 
(1979). For axisymm etric modes (k.<j> = 0), we re cover the 
dispersion relation of iGoldreich k, Schubertl |l967i ) without 
viscosity. 

We now set \ — and turn to evaluate the largest 
real part of the roots of the dispersion relation. We will 
seek the first order expansion of the growth rate in the form 
s = Nsq + flsi — N(sq + esi). For the largest real root at 
zeroth order, we get 



Close to marginal stability, the marginal root s m = 
without rotation can have the largest real part for slow rota- 
tion. However, s m is first order in e and the associated fluxes 
are of order e 2 and we safely neglect it. 

Finally, since so and si are always real, we simply take 

crfc = N(s + esi) (53) 
when iV 2 < and ah — otherwise. 

3.4 Convective fluxes 

Both the first and second order of the growth rate are real 
numbers, consequently the linear system of equations (|37|l to 
(|47|l introduces no phase shift between the perturbed fields 
involved. In our notations, the phase shift which enters the 
expression (|35p for the flux is ip(k) — or n for all pairs of 
variables of interest. 



3.4-1 Kinetic energy 

We start by deriving a useful relation between variables v' r 
and v'e- We use equation (|37p to express the variable v'^ in 
terms of the other two components of the velocity. Then, we 
combine equations (|39|l and (|40[) to eliminate the variable p' 
and so find the relationship between v' r and v'e to be 



= [(1 



[a/3s-l3-y2nM r + cry2tt cos 0]v' r 
a 2 )s-f3~f2£lMg + p-y2Qcosd]v'e. 



where we defined the more compact variables 

a =k.r = k x /k, 

P =k.O = ky/k and 

7 =k.(f> — k z /k. 

Here a 2 + /3 2 + 7 2 = 1. 

We now develop to first order in e the ratio of v'e over 
v' r from relation (|54|l to formally obtain 



(54) 

(55) 
(56) 
(57) 



v'e/v'r 



,2Tl2 + )' 



(58) 



/3 2 +7 2 

where / is a complicated function of a 2 and /3 2 of order 1. 
The terms involving si cancel. We now use equation (f37)l to 
get 



«7 



+ 7 e/(« ,n 



(59) 



so = 



(fe.f) 2 



(51) 



provided iV 2 < 0, which is the condition for instability. This 
expression shows that the fastest growing modes have zero 
radial wave number so that vertical convective plumes are 



/3 2 +7 2 

and we apply our saturation prescription (|27p to arrive at 

vl>k~ 2 ~k- 2n = \v'r\ 2 +\v'e\ 2 +\v'<i>\ 2 = v'l (l-4 e a/3 7 /)/(/3 2 +7 2 ). 

(60) 
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In the expression ()35[) we separate the integral over the mag- 
nitude of the wave vector from the integral over all possible 
directions of the wave vector. We determine that 



-r^-k 2n 4nk 2 dk 



2, a 2 . 2x dad/3d7 
(61) 



to the lowest (zeroth) order in e and we use k = k/k m . The 
first order in e is odd in 7 and its integration over the unit 
sphere yields a zero contribution. We perform the integral 
on the non-dimensional modulus of the wave vector k and 
write 



{v 2 } = 



47V 2 Xj 
Ti(2n - 3) 



with 



F rr = (4) s ~ 0.533 



(62) 



(63) 



where ()g denotes averaging over the unit sphere and we 



used So = 1 
(n = 11/6), 



/3 2 + 7 2 . 



With Kolomogorov scaling 



(v?) = —A 

7T 



2 N 2 F r , 



(64) 



Note that with Bolgiano-Obukhov scaling (n = 21/10), 
the prefactor 6/n decreases to 10/37T which is about twice 
smaller. We retain Kolomogorov scaling in the following. 
The other diagonal components of the Reynolds-stress ten- 
sor are 

6,2 .,2 „ 



of their value of 0.186. Our saturation prescription prob- 
ably overestimates the anisotropy because it neglects the 
tendency to isotropy at small scales down the turbulent cas- 
cade fcf. iRinconl l2006l ). Nevertheless, our model correctly 
accounts for the fact that the anisotropy does not depend 
on the latitude. 



3.4.2 Thermal fluxes 
From equation (|47|) . we write 



with 



and 



T' , 6 XjN 3 
1 ir g 



F Tr = (s 3 )s -0.589 



F T b = -{s a/3) s = 



(69) 
(70) 

(71) 



to the lowest order in e. The latitudinal thermal flux is zero 
to first order, consistent with the direction of the thermal 
gradients being vertical. 

The average equation for the evolution of thermal en- 
ergy is 



(65) 



Td t S + 'Ldr^Tr) - V.(pc pX VT) 
with the thermal convective flux given by 



(72) 



with 



and 



with 



F ee = (a 2 p 2 ) B 



{«/) = -A^iV 2 ^, 



= (a 7 >s = Fee ^ 0.067. 



(66) 
(67) 
(68) 



Our model predicts a strong anisotropic distribution of ve- 
locities with motions mostly in the radial direction. In accor- 
dance with the symetry of the problem, our model predicts 
cquipartition between the azimuthal and latitudinal direc- 
tions^ 

iKapyla et al.l |2004h compute these quantities in a num- 
ber of simulations of convection including rotation. Our 
small parameter e — Q/N translates in their notations as 
e = I ( P fla" ) 1 ■ Their simulation with Co — 1 corresponds 
to our e — 0.09. On the other hand we need the thermal dif- 
fusion timescale to be small before the rotation timescale be- 
cause we neglected thermal diffusion (and viscosity), and we 
require x/Hp/Vt to be small where H p is the pressure scale 
height. The value of this parameter is 0.17 for their simula- 
tion with Co = 1 and bigger for lower rotation rates, so we 
consider only their results at Co — 1. Using A m = \H P in 

equations {64]), (J65J) and (|6Tj) we find {v' 2 )? = 0.120 which 
is only slightly bigger than the value 0.090 which they find 
(viscous damping or a smaller A m could bring these values 
closer to one another). We also predict the ratio of horizon- 
tal to vertical motions {v'gVg + v'^v 1 ^) / (v' r v' r ) — 0.125 instead 



1 f* T' 
J r r = -pc p J dO sm6{v' r —}. 

Using equation (169[) we write 

_ 6 XlN 3 

J- r — pCp — tTr , 

TV q 



(73) 



(74) 



which we further develop into the more familiar thermal 
diffusive flux 



T r = — Dt pc p Ad r S, 



with the effective diffusion coefficient 



D T = - Ftt 

TV 



0.68NK 



(75) 



(76) 



MLTs traditionally make use of a diffusion coefficient of the 
form 



-Ok 



-Nt 2 > 

^ 1 v *rnj 



(77) 



where £ m ix is the mixing length. Comparing this expressions 
to the usual ID MLT, we can readily identify our smoothing 
length A m with the mixing length to a numerical factor of 

order one. 

In the simulations of IKapyla" et al l l|2004h with Co = 
1, they compute the eddy heat conductivity, Xrr = 
(v' r T' /T)g/N 2 in our notations. They compute the ratio 
Xrr/vt where v t = {v' 2 )^d/3 and d is the size of the con- 
vective zone (see their figure 19) and find it is between 0.5 
and 0.6 depending on the latitude. With A m = \H P , we pre- 
dict a slightly bigger value of 0.67 for this number, which is 
overestimated by about the same factor than for the r.m.s. 
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velocity (a lower value for A m would fix both numbers at the 
same time). 



3.4-3 Momentum fluxes 

The expressions for momentum fluxes are a bit more com - 
plicated. It has become common practice (sec Rudigcr 1989) 
to separate the momentum fluxes into a term linear in the 
rotation frequency (the A-effect) and a term which depends 
on the gradients of the rotation frequency (a-effect), rather 
than the gradients of specific angular momentum which we 
have used here. Therefore we offset the quantities M r and 
Me by their respective values for solid body rotation to com- 
pare more directly with previous work. Each momentum 
flux develops into a linear combination of terms, character- 
ized by four constant coefficients. For instance, we develop 
the radial momentum flux as 



I A3, Ml 



0.4 
0.2 
0.0 
-0.2 
-0.4 
-0.6 
-0.8 
-1.0 
-1.2 
-1.4 



1 1 1 1 


• 

• 






• • Kapyla (2004) 
- - KR93 

Gough (1978) 

Gough (2012} 
present work 





100 



F r <t, , M r ( M r ~ Sm 8 ) + F r 4, , MB ( Mo — COS 8 ) + F r 4, , cos 6 cos ( 



120 140 

Colatitude 



160 



180 



Since F, 



r<f>,M6 



= F„, 



— 0, this reduces to 



(v' r vL) = — A^iVfi (^Fr^Mr^rT— + F r d>, sin a) SIT18 

7r '2 omr 



with 

and 

-^7-0, sin 6 



-2 S g(a 2 + /3 2 )) s 



-0.687, 



(78) 



(79) 



(80) 



-2s /Hs 



-*o>s 



-0.589. (81) 



The two terms in expression (1791) represent respectively 
the a-effect and the A-effect ( see iRiidigerll 1989 ) . The radial 
a-effect is linked to differential ro tation and is diffusive in 
character. This was also found bv iHoudek fc Gougbl (|200lh 
at the equator (8 — 7r/2), but some quantities in their ex- 
pression are defined only implicitly which makes a direct 
comparison difficult. In the case of solid body rotation, we 
predict a A-effect in the form 



6 2 
A r< £ = — F rr f, tS i n e X m NQ sin t 

7T 



-l.l2\l,NQsm8. (82) 



This compares very well wit h the work of b oth 
IKichatinov fe Rudigerl (|l993T l and iGaraud et all (|2010h in 
the slow ro tation limit (see in particular equations (75) 
and (76) of IGaraud et al.ll2010l ) : How ever, we note the A- 
effect of IKichatinov fc RudigeiT (|l993T ) is essentially due to 
density gradients which we have neglected here, and they 
find a A-effect with the o pposite sign compared to us. Note 
that IGaraud et alj l|201Ch present their results as a function 
of anisotropy but, as they point out, the anisotropy is not 
arbitrary in their framework as in ours. In effect, the numer- 
ical coefficient in front of their A-term depends on the values 
of their closure parameters Ci, C2, C% and C7 and their ex- 
pressi on does not differ from that of IKichatinov fc Rudigerl 
(119931 ), except possibly for the numerical value of the pre- 
factor and the sign which could be either positive or nega- 
tive. Alth ough th e notio n of a A-term was not used a t that 
time, both IGoughl (|l978h and lDurnev fc Spruitl (|l979l ) have 
such a term in their formulation and correctly estimate its 

form in the slow ro tation limit. 

Simulations of lChanl (|200ll ), iRieutord et all (119941 ) and 



Figure 1. Normalised radial angular momentum flux due 
to the A-effect, (vIv'^ a/ ( v ' r v ' r )/ e i according to simulations 



in iKap yla et all (120041) (bl ue dots), to the predictions of 
iKichatinoT'fc'Ftudiged 1119931) (gree n) , to both predictions by 
iGoughl dl978l) (red, first order in e). lGoughl ||2012|) (cyan) and to 
our predictions (magenta). The numbers from the n umerical sim- 
ulatio ns are corrected from the shear (see table 3 in l Kapyla et al.1 
2004). To plot the results bv lKichatinov fc Rudigerl (11993^ . we use 
the effective viscosity vt defined as v± = (v' 2 )2 d/3 where d is the 
size of the convective region in their numeri cal setup and we too k 
(v' r v' r ) as measured in the simulations by Kapyla ct al. (2004). 
To plot the results by Gough, we used the parameter <J> = 1.9 as 
measured from the simulations. 



IKapyla" et al.l (|2004l ) all find a negative A-effect for slow ro- 
tation, in agreement with our results. However, we over- 
estimate by a large amount (up to a factor 4 near the equa- 
tor) the value of the transport co efficient compared to the 
simulations of IKapyla et al.1 12004), as seen in figure [1] The 
numbers extracted from the simulations are corrected from 
large scale shear flows which appear in their simulations. We 
use the numbers from their table 3, which shows the correc- 
tions themselves are of the same order as the measured radial 
momentum fluxes. 

In a similar way we write the latitudinal momentum 
flux as 



{veVi,) = -A m iVS2 -Fg^Me aa sin( 

TV Z Ou 



with 



-Ffl0,M9 = {—3OL P (7 



2 a 2 

- a p 



-0.123. 



(83) 



(84) 



The latitudinal a-effect is also diffusive but with a 
diffusion coefficient about six times smaller than the ra- 
dial one. The latitudinal A-effect for solid body rotation 
and our vertical entropy gradient is abse nt to first or - 
der in Q. This is in ag r eement with all of IGoughl Jl978h . 
Durnev fc Spruitl ll 19791 ). IKichatinov fc Rudigerl 1 19931) and 



Garaud et al. I |2010l ~ though in the case of lDurnev fc Spruitl 
l|l979l ) the latitudinal-azimuthal balance of kinetic energy is 
needed to cancel this te r m. T his is also consistent with the 
results of iKapyl a, et al.l (|2004t ) who find this term is much 
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smaller than its radial counterpart in the limit of slow rota- 
tion. 

The average equatio n for the trans port of angular mo- 
mentum can be found in IPurnevI jl985h (equation (3)). To 
first order the meridional circulation is absent and with our 
notations the angular momentum transport equation may 
be written as 

pd t r 2 sin 2 6 fi + ^^-d r (r 3 TZ rrt) ) + -^—de(sm 2 OTZse) = 
r x sm o 

(85) 

We now consider the special case of spherical symmetry 
which is more useful for ID stellar evolution codes. For this 
purpose we take Q to be a function of r only, so Mg = 
cos 8 and the latitudinal transport of momentum vanishes. 
It is however customary to integrate equation (|85[1 over the 
angles £ J™ d9 sin# x ... in order to eliminate the dg term 
so that 



p-r 2 d t O, + ^-d r (r z TZ r<t> ) = 0, 



with 



/ 

Jo 



Tin/, — d9 sin 9p(v r v,f,) 



(86) 



(87) 



d t r 2 Sl-\ — 3_9 r I pr 2 —F T &,Mr 
pr z ' 



dr{r ft) + 1 



F, 



refr,Mr 



(89) 

which shows that specific angular momentum is diffused 
with a diffusion coefficient 



— F^urXLN ~ 0.19 D T . (90) 

7T 



The A-effect yields an advection term which can be com- 
bined with the specific angular momentum gradient to pro- 
vide 



dtrQ. 



pr 



pr 2 D t 2q r 2 Q,d r In 



(91) 

which, after we evaluate the exponent of r in the loga- 
rithm, predicts a steady rotational profile in Q oc r -1 ' 86 . 
Thus the A-effect offsets the constant specific angular mo- 
mentum profile (SI oc r~ 2 ) by only a small amount. 
This contrasts with most ID studies of stellar rotation 
which assume solid body rotation in convection zones (e.; 
Mevnet fc Maederl Il997 ; iHeger, Langer fc Wooslevl l200o' . 
Potter. Tout fc Eldridgel (|201lT > studied the effects of vary- 
ing the specific angular momentum distribution in ID stellar 
models and found that the change in the total angular mo- 
mentum and additional shear generated at the boundary 
between convective and radiative regions can have a signifi- 
cant effect on the evolution of a star. 



4 CONCLUSION 

Using a generalized mixing length prescription, we have de- 
rived a self-consistent set of equations for axisymmetric 2D 
stellar evolution which includes a description of convective 
transport of angular momentum and heat. In the appendix 
A we list the full set of equations required to model the evo- 
lution of 2D stellar interiors at first order in Q./N as well 
as their ID spherically averaged equivalents. 

The thermal and momentum fluxes in radial and latitu- 
dinal directions are linked to the properties of the most un- 
stable local linear modes. In this respect our work in essence 



which also reads 



TLr4> = —pX m NQ / dO sin 6 [M r F r ^,,Mr + sin6> (F r j>, S m6 — F, 



r4> 



(88) 

We put the last expression back in the average momentum 
equation to obtain 



follows the spirit that iGougbl l| 19781) pioneered to estimate 
the fluxes due to small scale turbulent motions. However, 
our approach uses the angular directional dependence of the 
convective linear growth rate and determines the orientation 
of the convective motions. Thus, our prescription uses only 
one parameter, the smoothing length A m , which is readily 
seen to correspond to the mixing length in the ID limit. 
We have also studied the dynamics of convective motions in 
the presence of an arbitrary rotation field, with radial and 
latitudinal shear, as well as a radial and latitudinal thermal 
stratification. We provide simplified expressions relevant for 
special cases which can readily be incorporated in stellar 
evolution codes when the rotation is slow to first order in 
Q/N. The second order immediately brings features such as 
meridional circulation, non radial effective gravity and ther- 
mal gradients, and all terms of the dispersion relation need 
to be retained. In the future, we hope to be able incorpo- 
rate these ingredients in our formalism as well as to include 
magnetic fields. 
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APPENDIX A: SUMMARY OF STELLAR 
EVOLUTION EQUATIONS FOR SLOW 
ROTATION 

We reproduce here equations for the evolution of spherically 
symmetric slowly rotating stellar interiors, valid at first or- 
der in the rotation rate: 



-d r p + 3 = 
P 



(Al) 
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Coefficient 




Expression 


Value 


Frr 




(4)s 


0.533 


Foe = 




<a 2 /3 2 ) s 


0.067 


F Tr 




(4)s 


0.589 




(- 


-2 S 3(a 2 +/3 2 )> s 


-0.687 


^7-0, sin f? -^Tr 




(--8>s 


-0.589 




"a 


3 2 (7 4_ a2/3 2_ | g4 ))s 


-0.123 



Table Al. Coefficients relevant to the various correlations in- 
volved in the ffuxes. 

and 

TdtS + ^d r (r 2 D T pc p Ad r S) - V.(pc pX VT) = q (A2) 
with 

D T = - F T r N\ 2 m (A3) 

TV 

where N, the absolute magnitude of the square root of 

N 2 = ^ d r S (A4) 

is the buoyancy frequency and A m is our only parameter. We 
suggest to take the smoothing length A m as a given fraction 
of the pressure scale height as is usually done for the mix- 
ing length. Our comparison with numerical simulations and 
classical MLT suggests A m = \H V might be a reasonable 
choice. 

Poisson's equation reduces to 

g(r) = %{ Ar'A-Kpr' 2 (A5) 
r Jo 

and the usual boundary conditions are employed. The an- 
gular momentum evolution follows the equation 

pdtr 2 sin 2 0Q, + ^^9 r (r 3 1Z r d>) + — — ; rSsfsin 2 OTZse) = 
r £ sin o 

(A6) 

with 

Hr6 = —pX^NQ (^F r 6,Mr^T— + Fr6,suis) sin 6 (A7) 

7r / o In r 

and 

H<j>e = ^pX^Ntt ^Fe^.Affl ^^qJ 1 sinfi 1 - ( A 8) 

We summarize in table IaTI the linear coefficients Fj needed 
to determine the convective fluxes. 

When the rotation rate is taken spherically symetric, 
we obtain 

cV 2 fi - -^d r {pr 2 0.19D T r 2 Qd r In [r 1 ' 86 !^] } = (A9) 
for the transport of angular momentum. 
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